Neural network-based narma-l2 multivariable control method

ABSTRACT

The present invention discloses a neural network-based NARMA-L2 multivariable control method, which includes: deriving a NARMA-L2 multivariable control law of a nonlinear discrete system; using neural networks to offline identify nonlinear functions in the control law; designing the controller with the control law, designing a closed-loop control system for a certain type of turbofan engine, and adding the function of making the online update of neural network errors on the basis of this to correct the controller parameters, and studying the tracking performance of the control system. Aimed to minimize the quadratic performance index for the errors of the engine&#39;s high-pressure rotating speed and pressure ratio, by using neural networks in combination with the deduced NARMA-L2 multivariable control law, the present invention has designed a dual-variable controller, calculated the engine&#39;s control variables and controlled the fuel flow of the engine and the critical interface area of the exhaust nozzle.

TECHNICAL FIELD

The present invention belongs to the field of aero-engine simulation and control and specifically involves a neural network-based NARMA-L2 multivariable control method.

DESCRIPTION OF RELATED ART

The performance of early aero-engines was not high. They generally only needed to control the fuel flow to keep the engine speed constant. However, with higher requirements on the performance of aero-engines, the number of variables needed to control gradually increases in order to meet the control requirements of the aero-engines. For example, for a turbojet engine with afterburner, the afterburner fuel supply must also be used as a control variable in order to ensure that the nozzle pressure ratio of the turbine remains unchanged. It has been inevitable for the control system of the aero-engine to develop from a single-variable control system to a multi-variable control system.

There are limitations in the high-precision real-time modeling of aero-engines. Excessively accurate descriptions are not appropriate in the actual application of controllers. People often use linearized models in a small neighborhood of the equilibrium state to design controllers. This method satisfies the superposition principle, has a simple calculation process and can generally obtain satisfactory accuracy within the range of normal use. The nonlinear autoregressive-moving average-L2 (NARMA-L2) model has been applied in the field of engine control because it can accurately describe the input-output relationship of the nonlinear system, has a simple structure and can easily calculate the control law.

However, the currently commonly used NARMA-L2 control method is used for single-variable control. In the multivariable control process, decentralized loops are commonly used to form a multivariable control system, which is essentially just a combination of multiple single-variable control systems. This combined control method does not consider the impact of coupling and produces some impact on the control accuracy of the system.

SUMMARY

The technical problem to be solved by the present invention is the defect of the background technology. To solve this problem, the present invention has deduced the NARMA-L2 equation design control law of the multi-input multi-output (MIMO) system and designed a dual-variable controller involving the high-pressure rotating speed and pressure ratio of the turbofan engine to solve the coupling problem brought by the decentralized loop control system. Simulation results show that the controller has good stable and dynamic state performance.

The present invention adopts the following technical solution to solve the above technical problem:

A neural network-based NARMA-L2 multivariable control method, which includes the following steps:

Step A) Deduce the NARMA-L2 equations and multivariable control law of the MIMO nonlinear discrete system;

Step B) Design a dual-variable controller involving the high-pressure rotating speed and pressure ratio of a turbofan engine and perform performance verification at design points and non-design points and within the full flight envelope.

As a further optimization scheme for the neural network-based NARMA-L2 multivariable control method of the present invention, the specific sub-steps of Step A) are as follows:

Step A1) Obtain the NARMA equations of the MIMO system through recursion according to the state space description of the MIMO nonlinear discrete system;

Step A2) According to the NARMA equations of the MIMO system, perform multivariate Taylor expansion near the equilibrium point and ignore the Taylor high-order remainders above the quadratic term to obtain the NARMA-L2 equations of the MIMO system;

Step A3) Operate the NARMA-L2 equations of the dual-input dual-output system to obtain the dual-variable control law of the system;

As a further optimization scheme for the neural network-based NARMA-L2 multivariable control method of the present invention, the specific sub-steps of Step B) are as follows:

Step B1) Take the design of the dual-variable controller of a certain type of turbofan engine as an example. Six non-linear functions in the dual-variable control law are identified offline through six neural networks. The network model adopts the BP-NN network topology and uses the input and output data at the past moment as the input of the neural network and the mapping result of the nonlinear function as the output of the neural network to complete the training;

Step B2) According to the control law, construct a matrix with the outputs of the nonlinear functions identified by the neural network, and establish a closed-loop controller for the high-pressure rotating speed and pressure ratio of the engine, get the values of the control variables online, and set up a quadratic performance index according to the error between the actual output and the command output to optimize neural network parameters online;

Step B3) After designing the dual-variable controller, establish the dual-variable closed-loop control system of the turbofan engine and verify its tracking performance at design points and non-design points and within the full flight envelope.

The neural network-based NARMA-L2 dual-variable controller of turbofan engine includes a neural network module and an online optimization module.

The above-mentioned neural network module is used to approximate the nonlinear functions in the control law according to the past input and output data of the system.

The above-mentioned online optimization module is used for the online correction of controller parameters. After the dual-variable controller obtains the controlled variable command at each moment, it calculates the control variable corresponding to the command according to the control law and then inputs it into the engine. The engine then outputs the real value of the controlled variable. In combination with the quadratic performance index, the controlled variable command and the actual output are used for the online correction of the neural network that approximates the nonlinear functions so as to track the trajectory of the engine within the full flight envelope.

Compared with the existing technologies, the above technical solution adopted by the present invention has the following technical effects:

(1) The NARMA-L2 multi-variable control method proposed by the present invention can form a single-loop multi-variable closed-loop control system and can eliminate the coupling impact of multiple variables and improve the control quality;

(2) The neural network-based NARMA-L2 dual-variable controller proposed by the present invention can track any command stably in the full flight envelope, and has good stable and dynamic state performance and shortens the response time of the controller. This verifies the effectiveness of the control method.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is the offline identification structure diagram of NARMA-L2 neural network;

FIG. 2 is the off-line identification error curve diagram of the neural network;

FIG. 3 is the structure diagram of NARMA-L2 dual-variable closed-loop control system;

FIG. 4 is the tracking response diagram of the rotating speed outputted at the design point (H=0 m and Ma=0);

FIG. 5 is the tracking response diagram of the pressure ratio outputted at the design point (H=0 m and Ma=0);

FIG. 6 is the tracking response diagram of the rotating speed outputted at the non-design point (H=4,000 m and Ma=0.8);

FIG. 7 is the tracking response diagram of the pressure ratio outputted at the non-design point (H=4,000 m and Ma=0.8);

FIG. 8 is the change diagram of altitude and Mach number during the take-off, climb, descent, acceleration and deceleration processes of the aircraft;

FIG. 9 is the tracking response diagram of the rotating speed outputted during the take-off, climb, descent, acceleration and deceleration processes of the aircraft;

FIG. 10 is the tracking response diagram of the pressure ratio outputted during the take-off, climb, descent, acceleration and deceleration processes of the aircraft.

DESCRIPTION OF THE EMBODIMENTS

The technical solution of the present invention will be further illustrated in detail below in conjunction with the accompanying drawings.

The idea of the present invention is first to obtain the NARMA equations of the MIMO nonlinear discrete system through recursion and to perform the multivariate Taylor expansion at the equilibrium point to obtain the NARMA-L2 equations of the system, and then to get the control law from the equations through simple derivation. In response to the multi-variable control requirements of the turbofan engine, the dual-variable controller is designed by using the derived NARMA-L2 multi-variable control law to approximate the nonlinear functions in the control law based on the neural network on the basis of the component-level model of the turbofan engine. Compared with the decentralized loop dual-variable control system, the controller eliminates the impact of coupling and has better tracking effects and a little higher accuracy. Compared with other multivariable control methods, this control method has a simple structure and can easily design the control law, and also has good stable and dynamic state performance.

The embodiment of the present invention takes the design of a dual-variable controller involving the high-pressure rotating speed and pressure ratio of a certain type of turbofan engine as an example. The NARMA-L2 multivariable control method includes the following steps:

Step A) Deduce the NARMA-L2 equations and multivariable control law of the MIMO nonlinear discrete system;

Step B) Design a NARMA-L2 dual-variable controller and perform performance verification at design points and non-design points and within the full flight envelope.

Wherein, the specific sub-steps of Step A) are as follows:

Step A1) The state space of the MIMO nonlinear discrete system is described as:

$\begin{matrix} {{x\left\lbrack {k + 1} \right\rbrack} = {f\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack}} \right)}} \\ {{y_{1}\lbrack k\rbrack} = {h_{1}\left( {x\lbrack k\rbrack} \right)}} \\ {{y_{2}\lbrack k\rbrack} = {h_{2}\left( {x\lbrack k\rbrack} \right)}} \\ \ldots \\ {{y_{m}\lbrack k\rbrack} = {h_{m}\left( {x\lbrack k\rbrack} \right)}} \end{matrix}$

Wherein, the input variable u[k]=[u₁[k]u₂[k]. . . u_(m)[k]]ϵR^(m), the state vector x[k]ϵR^(m), the output values, y₁[k]ϵR, y₂[k]ϵR, . . . , y_(m)[k]ϵR, the functions f(.),h₁(.), . . . , h_(m)(.)ϵC^(∞), and the original point is the equilibrium state, and k is the time.

In practical applications, because usually only input and output signals are obtained, it is needed to establish the input and output description of the system. Get the following results through recursion:

$\begin{matrix} {{y_{1}\lbrack k\rbrack} = {{h_{1}\left( {x\lbrack k\rbrack} \right)} = {\Psi_{1}\left( {x\lbrack k\rbrack} \right)}}} \\ {{y_{1}\left\lbrack {k + 1} \right\rbrack} = {{h_{1}\left( {x\left\lbrack {k + 1} \right\rbrack} \right)} = {{h_{1}\left( {f\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack}} \right)} \right)} = {\Psi_{2}\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack}} \right)}}}} \\ {{y_{1}\left\lbrack {k + 2} \right\rbrack} = {{h_{1}\left( {x\left\lbrack {k + 2} \right\rbrack} \right)} = {h_{1}\left( {f\left( {{x\left\lbrack {k + 1} \right\rbrack},{u\left\lbrack {k + 1} \right\rbrack}} \right)} \right)}}} \\ {\left. {= {h_{1}\left( {{f\left( {{f\left( {x\lbrack k\rbrack} \right)},{u\lbrack k\rbrack}} \right)},{u\left\lbrack {k + 1} \right\rbrack}} \right)}} \right) = {\Psi_{3}\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack}} \right)}} \\ \ldots \\ {{y_{1}\left\lbrack {k + n - 1} \right\rbrack} = {{h_{1} \circ {f^{n - 1}\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack}} \right)}} = {\Psi_{n}\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack},\ldots,{u\left\lbrack {k + n - 2} \right\rbrack}} \right)}}} \end{matrix}$

Define:

y_(1n)[k]=[y₁[k], y₁[k+1], . . . , y₁[k+n−1]]=Ψ(x[k], u_(n−1)[k]) u_(n−1)[k]=[u[k],u[k+1], . . . , u[k+n−2]]

Then the following equation can be deduced:

Ψ(x[k], u_(n−1)[k])=Y_(1n)[k]

Supposed that ∂Ψ/∂x is nonsingular, due to the existence theorem of the implicit function, it can be deduced that the state vector x[k] can be expressed with Y_(1n)[k] and u_(n−1)[k].

Also:

$\begin{matrix} {{x\left\lbrack {k + n} \right\rbrack} = {f\left( {{x\left\lbrack {k + n - 1} \right\rbrack},{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \\ {= {f\left( {{f\left( {{x\left\lbrack {k + n - 2} \right\rbrack},{u\left\lbrack {k + n - 2} \right\rbrack}} \right)},{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \\ {= \ldots} \\ {= {f \circ {f^{n - 1}\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack}} \right)}}} \\ {= {\Phi\left( {{x\lbrack k\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \\ {= {\Omega\left( {{Y_{1n}\lbrack k\rbrack},{u_{n - 1}\lbrack k\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \\ {= {\Omega\left( {{Y_{1n}\lbrack k\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \end{matrix}$

Wherein, Φ(.) and Ω(.) are smoothing functions.

According to y₁[k+n]=h₁(x[k+n]), the following can be obtained:

$\begin{matrix} {{y_{1}\left\lbrack {k + n} \right\rbrack} = {h_{1}\left( {x\left\lbrack {k + n} \right\rbrack} \right)}} \\ {= {h_{1}\left( {\Omega\left( {{y_{1}\lbrack k\rbrack},{y_{1}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{1}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)} \right)}} \\ {= {F\left( {{y_{1}\lbrack k\rbrack},{y_{1}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{1}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \end{matrix}$

Wherein, F(.), Φ(.) and Ω(.) are smoothing functions.

The following can be obtained by the same way:

$\begin{matrix} {{y_{2}\left\lbrack {k + n} \right\rbrack} = {h_{2}\left( {x\left\lbrack {k + n} \right\rbrack} \right)}} \\ {= {h_{2}\left( {\psi\left( {{y_{2}\lbrack k\rbrack},{y_{2}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{2}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)} \right)}} \\ {= {G\left( {{y_{2}\lbrack k\rbrack},{y_{2}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{2}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \end{matrix}$

$\begin{matrix} {{y_{m}\left\lbrack {k + n} \right\rbrack} = {h_{m}\left( {x\left\lbrack {k + n} \right\rbrack} \right)}} \\ {= {h_{m}\left( {P\left( {{y_{m}\lbrack k\rbrack},{y_{m}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{m}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)} \right)}} \\ {= {H\left( {{y_{m}\lbrack k\rbrack},{y_{m}\left\lbrack {k + 1} \right\rbrack},\ldots,{y_{m}\left\lbrack {k + n - 1} \right\rbrack},{u\lbrack k\rbrack},{u\left\lbrack {k + 1} \right\rbrack},\ldots,{u\left\lbrack {k + n - 1} \right\rbrack}} \right)}} \end{matrix}$

Wherein, ψ(.), G(.), P(.) and H(.) are smoothing functions.

Therefore, the NARMA equations of the MIMO system are recursively obtained as:

$\begin{matrix} {{y_{1}\left\lbrack {k + 1} \right\rbrack} = {F\left( {{y_{1}\left\lbrack {k - n + 1} \right\rbrack},{y_{1}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{1}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \\ {{y_{2}\left\lbrack {k + 1} \right\rbrack} = {G\left( {{y_{2}\left\lbrack {k - n + 1} \right\rbrack},{y_{2}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{2}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \\ \ldots \\ {{y_{m}\left\lbrack {k + 1} \right\rbrack} = {H\left( {{y_{m}\left\lbrack {k - n + 1} \right\rbrack},{y_{m}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{m}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \end{matrix}$

Wherein, the input variable u[k]=[u₁[k]u₂[k]. . . u_(m)[k]]ϵR^(m), the output values y₁[k]ϵR, y₂[k]ϵR, . . . , y_(m)[k]ϵR. the functions, F(.), G(.), H(.)ϵC^(∞) and k is the time.

Step A2) The NARMA model is actually an accurate mathematical expression of the nonlinear system in the equilibrium state, but it is not suitable for the actual system due to a large amount of calculation, so the NARMA model needs to be approximated and simplified.

Make:

M₁ = (y₁[k − n + 1], y₁[k − n + 2] , …, y₁[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1])M₂ = (y₂[k − n + 1], y₂[k − n + 2] , …, y₂[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1])…M_(m) = (y_(m)[k − n + 1], y_(m)[k − n + 2] , …, y_(m)[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1])

By performing multivariate Taylor expansion for the NARMA equations near the equilibrium point and ignoring Taylor high-order remainders, the NARMA-L2 equations of the system are obtained as:

$\begin{matrix} {{y_{1}\left\lbrack {k + 1} \right\rbrack} = {F\left( {M_{1},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {F\left( {M_{1},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial F}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial F}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial F}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \end{matrix}$ $\begin{matrix} {{y_{2}\left\lbrack {k + 1} \right\rbrack} = {G\left( {M_{2},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {G\left( {M_{2},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial G}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial G}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial G}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \\ \ldots \end{matrix}$ $\begin{matrix} {{y_{m}\left\lbrack {k + 1} \right\rbrack} = {G\left( {M_{m},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {G\left( {M_{m},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial G}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial G}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial G}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \end{matrix}$

Step A3) Convert the NARMA-L2 equations of the MIMO system into a matrix form:

$\begin{bmatrix} {y_{1}\left( {k + 1} \right)} \\ {y_{2}\left( {k + 1} \right)} \\  \vdots \\ {y_{m}\left( {k + 1} \right)} \end{bmatrix} = {\begin{bmatrix} F_{0} \\ G_{0} \\  \vdots \\ H_{0} \end{bmatrix} + {\begin{bmatrix} F_{1} & F_{2} & \ldots & F_{m} \\ G_{1} & G_{2} & \ldots & G_{m} \\  \vdots & \vdots & & \vdots \\ H_{1} & H_{2} & \ldots & H_{m} \end{bmatrix}\begin{bmatrix} {u_{1}(k)} \\ {u_{2}(k)} \\  \vdots \\ {u_{m}(k)} \end{bmatrix}}}$

Wherein,

${{{{{{{{{{{{{{{{{F_{0} = F},{F_{1} = \frac{\partial F}{\partial{u_{1}(k)}}}}❘}_{({M_{1},0,0,\ldots,0})},{{F_{2} = \frac{\partial F}{\partial{u_{2}(k)}}}❘}_{({M_{1},0,0,\ldots,0})},{F_{m} = \frac{\partial F}{\partial{u_{m}(k)}}}}❘}_{({M_{1},0,0,\ldots,0})},{G_{0} = G},{G_{1} = \frac{\partial G}{\partial{u_{1}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{G_{2} = \frac{\partial G}{\partial{u_{2}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{G_{m} = \frac{\partial G}{\partial{u_{m}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{H_{0} = H},{H_{1} = \frac{\partial H}{\partial{u_{1}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})},{H_{2} = \frac{\partial H}{\partial{u_{2}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})},{H_{m} = \frac{\partial H}{\partial{u_{m}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})}$

represent a nonlinear function respectively, and their input parameters are (y₁[k−n+1], y₁[k−n+2], . . . , y₁[k], u[k−1], . . . , u[k−n+1]).

Therefore, according to the matrix operation, the multivariable control law is obtained as:

$\begin{bmatrix} {u_{1}(k)} \\ {u_{2}(k)} \\  \vdots \\ {u_{m}(k)} \end{bmatrix} = {\begin{bmatrix} F_{1} & F_{2} & \ldots & F_{m} \\ G_{1} & G_{2} & \ldots & G_{m} \\  \vdots & \vdots & & \vdots \\ H_{1} & H_{2} & \ldots & H_{m} \end{bmatrix}^{- 1}\left( {\begin{bmatrix} {y_{1}\left( {k + 1} \right)} \\ {y_{2}\left( {k + 1} \right)} \\  \vdots \\ {y_{m}\left( {k + 1} \right)} \end{bmatrix} - \begin{bmatrix} F_{0} \\ G_{0} \\  \vdots \\ H_{0} \end{bmatrix}} \right)}$

The specific sub-steps of Step B) are as follows:

Step B1) Take the fuel flow of the turbofan engine and the critical cross-sectional area of the exhaust nozzle as control variables, and take the high-pressure rotating speed and pressure ratio of the engine as controlled variables to design a dual-variable controller. For the dual-variable control system, its control law is:

$\begin{bmatrix} {u_{1}(k)} \\ {u_{2}(k)} \end{bmatrix} = {\begin{bmatrix} F_{1} & F_{2} \\ G_{1} & G_{2} \end{bmatrix}^{- 1}\left( {\begin{bmatrix} {y_{1}\left( {k + 1} \right)} \\ {y_{2}\left( {k + 1} \right)} \end{bmatrix} - \begin{bmatrix} F_{0} \\ G_{0} \end{bmatrix}} \right)}$

Wherein, non-linear functions F0, F1, F2, G0, G1, G2 can be fitted with a neural network, respectively. The six nonlinear functions in the dual-variable control law are identified offline through six neural networks, and the network model adopts the BP-NN network topology.

The three neural networks approximating the nonlinear functions F0, F1 , F2 all take M₁=(y₁[k−n+1], y₁[k−n+2], . . . , y₁[k], u[k−n+1], u[k−n+2], . . . , u[k−1]) as the input, and their outputs are the approximate values of F₀(y₁[k−n+1], y₁[k−n+2], . . . , y₁[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), F₁(y₁[k−n+1], y₁[k−n+2], . . . , y₁[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), F₂(y₁[k−n+1], y₁[k−n+2], . . . , y₁[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), respectively. The turbofan engine is generally recognized as a second-order system, so n is equal to 3, the input neuron of the neural network is 9, the hidden layer's neuron is 25, the output layer's neuron is 1 and the hidden layer's activation function is the tanh function.

The three neural networks approximating the nonlinear functions G0, G1, G2 all take M₂=(y₂[k−n+1], y₂[k−n+2], . . . , y₂[k], u[k−n+1], u[k−n+2], . . . , u[k−1]) as the input, and their outputs are the approximate values of G₀(y₂[k−n+1], y₂[k−n+2], . . . , y₂[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), G₁(y₂[k−n+1], y₂[k−n+2], . . . , y₂[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), G₂=(y₂[k−n+1], y₂[k−n+2], . . . , y₂[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), respectively. Also, n is equal to 3, the input neuron of the neural network is 9, the hidden layer's neuron is 25, the output layer's neuron is 1 and the hidden layer's activation function is the tanh function.

The engine model is used to generate a large amount of ground point input and output data. Such data are used as the input of each neural network. The model output of the NARMA-L2 equation is calculated through the output of each neural network, and the error is backpropagated to train each neural network.

The performance index of the model is:

${J = {\frac{1}{2}\left( {{e_{1}^{2}\lbrack k\rbrack} + {e_{2}^{2}\lbrack k\rbrack}} \right)}}{e_{1} = {{y_{1}\lbrack k\rbrack} - {y{m_{1}\lbrack k\rbrack}}}}{e_{2} = {{y_{2}\lbrack k\rbrack} - {y{m_{2}\lbrack k\rbrack}}}}$

Wherein, ym₁ and ym₂ respectively represent the model outputs calculated according to the NARMA-L2 equations, and y₁ and y₂ respectively represent the actual high-pressure rotating speed and pressure ratio of the engine.

The offline identification structure of the NARMA-L2 model is shown in FIG. 1, where the bolded variables are all matrices, namely u [k]=[u₁[k]u₂[k]]^(T), y[k]=[y₁[k]y₂[k]]^(T), ym[k]=[ym[k]ym₂[k]]^(T), e=[e₁ e₂ ]^(T). The NARMA-L2 model identification error curve is shown in FIG. 2. As the number of training increases, the neural network identification is basically completed. The theoretical output values calculated with the NARMA-L2 equations are basically consistent with the actual output values.

Step B2) According to the NARMA-L2 dual-variable control law, construct a matrix with the outputs of the nonlinear functions already identified by the neural networks, establish a closed-loop controller for the high-pressure rotating speed and pressure ratio of the engine, and get the values of the control variables online. Moreover, establish a quadratic performance index according to the error between the actual output and the command output, and use the backpropagation algorithm to optimize the neural network parameters online.

The performance index is:

${J = {\frac{1}{2}\left( {{e_{c1}^{2}\lbrack k\rbrack} + {e_{c2}^{2}\lbrack k\rbrack}} \right)}}{e_{c1} = {{y_{1}\lbrack k\rbrack} - {r_{1}\lbrack k\rbrack}}}{e_{c2} = {{y_{2}\lbrack k\rbrack} - {r_{2}\lbrack k\rbrack}}}$

Wherein, r₁ and r₂ represent the engine's high-pressure rotating speed and pressure ratio commands, respectively; and y₁ and y₂ represent the engine's actual high-pressure rotating speed and pressure ratio, respectively.

The built control structure is shown in FIG. 3, where the bolded variables are all matrices, namely u[k]=[u₁[k]u₂[k]]^(T), y[k]=[y₁[k]y₂[k]]^(T), r[k]=[k]r₁[k]r₂[k]]^(T), e_(c)=[e _(c1) e_(c2)]^(T).

Step B3) After the dual-variable controller is designed, the dual-variable closed-loop control system of the turbofan engine is established and its tracking performance is verified at design points and non-design points and within the full flight envelope.

The simulation results under the state of the design point (H=0 m and Ma=0) are shown in FIG. 4 and FIG. 5. It can be seen from the figures that the high-pressure rotating speed and pressure ratio under this state can track the command step signal stably without steady-state error; the overshoot of the high-pressure rotating speed is 6.32% and the adjustment time is less than 2 s; the overshoot of the pressure ratio is 5.53% and the adjustment time is less than 1.4 s, so it has good stable and dynamic state performance.

The simulation results at the non-design point (H=4,000 m and Ma=0.8) are shown in FIG. 6 and FIG. 7. It can be seen from the figures that the high-pressure rotating speed and pressure ratio under this state can track the change of the command signal stably, basically without steady-state error; the overshoot of the high-pressure rotating speed is 7.95% and the adjustment time is less than 2.5 s; and the overshoot of the pressure ratio is 6.34% and the adjustment time is less than 2 s, so it has good dynamic characteristics.

The changes of the altitude and Mach number during the take-off, climb, descent, acceleration and deceleration process of the aircraft are shown in FIG. 8, and the output tracking effects of the high-pressure rotating speed and pressure ratio are shown in FIG. 9 and FIG. 10. It can be seen from the figures that in various processes, the high-pressure rotating speed and pressure ratio can basically track the change of the command stably, and the system responds quickly, so the model can well meet the needs of the dual-variable control, has good control quality and can adapt to the changes of various working conditions of the engine.

It should be pointed out that what mentioned above is only the embodiment of the present invention, but the protection scope of the present invention is not limited to this. Any changes and replacements that any technician familiar with this technical field can easily make within the technical scope disclosed by the present invention should fall into the protection scope of the present invention. Therefore, the protection scope of the present invention shall be the protection scope of the following claims. 

1. A neural network-based NARMA-L2 multivariable control method, comprising following steps: Step A) deducing the NARMA-L2 equations and multivariable control law of a multi-input multi-output (MIMO) nonlinear discrete system; Step B) using neural networks to offline identify nonlinear functions in the muitivariable control law; using the identified neural networks to design a dual-variable controller, establishing a dual-variable closed-loop control system of a turbofan engine, and verifying a stable and dynamic state performance of the dual-variable controller, wherein specific sub-steps of the Step A) are as follows: Step A1) obtaining NARMA equations of the MIMO system through recursion according to a state space description of a nonlinear discrete system, and apply it to the MIMO system; Step A2) according to the NARMA equations of the MIMO system, performing multivariate Taylor expansion near an equilibrium point and ignore Taylor high-order remainders to obtain the NARMA-L2 equations of the NEMO nonlinear system; Step A3) performing matrix operations by using the NARMA-L2 equations of the MIMO nonlinear system to obtain the multivariable control law of the MIMO system, wherein the state space description of the MIMO nonlinear discrete system in the Step A1) is as follows: x[k + 1] = f(x[k], u[k])y₁[k] = h₁(x[k])y₂[k] = h₂(x[k])⋯y_(m)[k] = h_(m)(x[k]) wherein, u[k]=[u₁[k]u₂[k]. . . . u_(m)[k]]ϵR^(m) is a control variable input of the system, y₁[k]ϵR, y₂[k]ϵR, . . . , y_(m)[k]ϵR is a controlled variable output of the system, m represents a number of the selected control variables and controlled variables, a state vector x[k]ϵR^(n), n represents a dimension of the state vector, functions f(.), h₁(.), . . . , h_(m)(.) ϵC^(∞), and an original point is an equilibrium state, and k represents time; the NARMA equations of the MIMO system are: $\begin{matrix} {{y_{1}\left\lbrack {k + 1} \right\rbrack} = {F\left( {{y_{1}\left\lbrack {k - n + 1} \right\rbrack},{y_{1}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{1}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \\ {{y_{2}\left\lbrack {k + 1} \right\rbrack} = {G\left( {{y_{2}\left\lbrack {k - n + 1} \right\rbrack},{y_{2}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{2}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \\ \ldots \\ {{y_{m}\left\lbrack {k + 1} \right\rbrack} = {H\left( {{y_{m}\left\lbrack {k - n + 1} \right\rbrack},{y_{m}\left\lbrack {k - n + 2} \right\rbrack},\ldots,{y_{m}\lbrack k\rbrack},{u\left\lbrack {k - n + 1} \right\rbrack},{u\left\lbrack {k - n + 2} \right\rbrack},\ldots,{u\lbrack k\rbrack}} \right)}} \end{matrix}$ wherein, an input variable a u[k]=[u₁[k]u₂[k]. . . . u_(m)[k]]ϵR^(m), output values y₁[k]ϵR, y₂[k]ϵR, . . . , y_(m)[k]ϵR, functions F(.), G(.), H(.)ϵC^(∞) and k represents time, wherein the NARMA-L2 equations of the MIMO system in the Step A2) are: $\begin{matrix} {{y_{1}\left\lbrack {k + 1} \right\rbrack} = {F\left( {M_{1},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {F\left( {M_{1},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial F}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial F}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial F}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{1},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \end{matrix}$ $\begin{matrix} {{y_{2}\left\lbrack {k + 1} \right\rbrack} = {G\left( {M_{2},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {G\left( {M_{2},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial G}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial G}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial G}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{2},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \\ \ldots \end{matrix}$ $\begin{matrix} {{y_{m}\left\lbrack {k + 1} \right\rbrack} = {H\left( {M_{2},{u_{1}\lbrack k\rbrack},{u_{2}\lbrack k\rbrack},\ldots,{u_{m}\lbrack k\rbrack}} \right)}} \\ {= {H\left( {M_{m},0,0,\ldots,0} \right)}} \\ {{{{{{{{{+ \frac{\partial H}{\partial{u_{1}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{1}\lbrack k\rbrack}} + \frac{\partial H}{\partial{u_{2}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{2}\lbrack k\rbrack}} + \ldots + \frac{\partial H}{\partial{u_{m}\lbrack k\rbrack}}}❘}_{({M_{m},0,0,\ldots,0})}{u_{m}\lbrack k\rbrack}} \end{matrix}$ wherein, M₁ = (y₁[k − n + 1], y₁[k − n + 2] , …, y₁[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1])M₂ = (y₂[k − n + 1], y₂[k − n + 2] , …, y₂[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1])…M_(m) = (y_(m)[k − n + 1], y_(m)[k − n + 2] , …, y_(m)[k], u[k − n + 1], u[k − n + 2] , …, u[k − 1]).
 2. (canceled)
 3. The neural network-based NARMA-L2 multivariable control method according to claim 1, wherein specific sub-steps of the Step B) are as follows: Step B1) first, taking the fuel flow of the turbofan engine and a critical cross-sectional area of an exhaust nozzle as control variables, and taking a high-pressure rotating speed and a pressure ratio of the turbofan engine as controlled variables, using an engine model to generate input and output data, and using neural networks to offline identify the nonlinear functions in the control law; Step B2) according to the control law, constructing a matrix with outputs of the nonlinear functions already identified by the neural networks, establishing a closed-loop controller for the high-pressure rotating speed and the pressure ratio of the turbofan engine, and obtain values of the control variables online, establishing quadratic performance index according to an error between actual output and a command output, and optimizing neural network parameters online; Step B3) after the dual-variable controller is designed, verifying a tracking performance of the dual-variable controller at design points and non-design points and within a full flight envelope. 4-5. (canceled)
 6. The neural network-based NARMA-L2 multivariable control method according to claim 1, wherein the multivariable control law of the system in the Step A3) is: $\begin{bmatrix} {u_{1}(k)} \\ {u_{2}(k)} \\  \vdots \\ {u_{m}(k)} \end{bmatrix} = {\begin{bmatrix} F_{1} & F_{2} & \ldots & F_{m} \\ G_{1} & G_{2} & \ldots & G_{m} \\  \vdots & \vdots & & \vdots \\ H_{1} & H_{2} & \ldots & H_{m} \end{bmatrix}^{- 1}\left( {\begin{bmatrix} {y_{1}\left( {k + 1} \right)} \\ {y_{2}\left( {k + 1} \right)} \\  \vdots \\ {y_{m}\left( {k + 1} \right)} \end{bmatrix} - \begin{bmatrix} F_{0} \\ G_{0} \\  \vdots \\ H_{0} \end{bmatrix}} \right)}$ ${{{{{{{{{{{{{{F_{m} = \frac{\partial F}{\partial{u_{m}(k)}}}❘}_{({M_{1},0,0,\ldots,0})},{G_{0} = G},{G_{1} = \frac{\partial G}{\partial{u_{1}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{G_{2} = \frac{\partial G}{\partial{u_{2}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{G_{m} = \frac{\partial G}{\partial{u_{m}(k)}}}}❘}_{({M_{2},0,0,\ldots,0})},{H_{0} = H},{H_{1} = \frac{\partial H}{\partial{u_{1}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})},{H_{2} = \frac{\partial H}{\partial{u_{2}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})},{H_{m} = \frac{\partial H}{\partial{u_{m}(k)}}}}❘}_{({M_{m},0,0,\ldots,0})}.$
 7. The neural network-based NARMA-L2 multivariable control method according to claim 3, wherein a function of the quadratic performance index established in the Step B2) is: ${J = {\frac{1}{2}\left( {{e_{c1}^{2}\lbrack k\rbrack} + {e_{c2}^{2}\lbrack k\rbrack}} \right)}}{e_{c1} = {{y_{1}\lbrack k\rbrack} - {r_{1}\lbrack k\rbrack}}}{e_{c2} = {{y_{2}\lbrack k\rbrack} - {r_{2}\lbrack k\rbrack}}}$ wherein, r₁ and r₂ represent high-pressure rotating speed and commands of the pressure ratio of the turbofan engine, respectively; and y₁ and y₂ represent actual high-pressure rotating speed and actual pressure ratio of the turbofan engine, respectively; after the dual-variable controller obtaining a command of controlled variable at each moment, the dual-variable controller calculates the control variable corresponding to the command according to the control law and inputs the control variable into the engine, the engine then outputs a real value of the controlled variable, in combination with the quadratic performance index, the command of the controlled variable and an actual output value are used for an online correction of the neural network that approximates the nonlinear function, so as to achieve the adaptive control of the engine in the full flight envelope. 